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ABSTRACT 

Context. The W51 complex hosts the supernova remnant W51C which is known to interact with the molecular clouds in the star forming region 
W51B. In addition, a possible pulsar wind nebula CXO J1923 18.5+ 140305 was found likely associated with the supernova remnant. Gamma-ray 
emission from this region was discovered by FermifLAT (between 0.2 and 50 GeV) and H.E.S.S. (>1 TeV). The spatial distribution of the events 
could not be used to pinpoint the location of the emission among the pulsar wind nebula, the supernova remnant shell and/or the molecular cloud. 
However, the modeling of the spectral energy distribution presented by the FermifLAT collaboration suggests a hadronic emission mechanism. 
The possibility that the gamma-ray emission from such an object is of hadronic origin can contribute to solvingthe long-standing problem of the 
contribution to galactic cosmic rays by supernova remnants. 

Aims. Our aim is to determine the morphology of the very-high-energy gamma-ray emission of W5 1 and measure its spectral properties. 
Methods. We performed observations of the W5 1 complex with the MAGIC telescopes for more than 50 hours. The energy range accessible with 
MAGIC extends from 50 GeV to several TeV, allowing for the first spectral measurement at these energies. In addition, the good angular resolution 
in the medium (few hundred GeV) to high (above 1 TeV) energies allow us to perform morphological studies. We look for underlying structures 
by means of detailed morphological studies. Multi-wavelength data from this source have been sampled to model the emission with both leptonic 
and hadronic processes. 

Results. We detect an extended emission of very-high-energy gamma rays, with a significance of 1 1 standard deviations. We extend the spectrum 
from the highest FermifLAT energies to ~ 5 TeV and find that it follows a single power law with an index of 2.58 ± 0.07 s t a t ± 0.22 svs t. The main 
part of the emission coincides with the shocked cloud region, while we find a feature extending towards the pulsar wind nebula. The possible 
contribution of the pulsar wind nebula, assuming a point-like source, shows no dependence on energy and it is about 20% of the overall emission. 
The broad band spectral energy distribution can be explained with a hadronic model that implies proton acceleration above 100 TeV. This result, 
together with the morphology of the source, tentatively suggests that we observe ongoing acceleration of ions in the interaction zone between 
supernova remnant and cloud. 

Key words. Acceleration of particles - cosmic rays - ISM: supernova remnants - ISM: clouds - Gamma rays: general - Gamma rays: ISM 
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1. Introduction 

W51 is a massive molecular complex located at the tangen- 
tial point (I = 49°) of the Sagitarius arm of the Galaxy, at 
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a distance of ~ 5.5kpc (ISato et alJ 12010). As seen in radio 
continuum images, three main components are identified: the 
star-forming regions W51A and W51B and, attached to the 
south-eastern boundary of W5 IB, the supernova remn ant (SNR) 
W51C. The estimated age of this SNR is 30 kyrs dKoo et al l 
1995a). Evidence of interaction between W51C and W51B is 
provided by several observations. Most cruc ial of them are the 
existence of two 1720 MHz OH masers dGreen etaLI 1 1997b 
and the detection of about 10 3 solar masses of atomic gas at 
a velocity shifted bet ween 20 and 120 km s" 1 with respect to 
its ambient medium dKoo & Moonl |1997ab . The high- velocity 
atomic ga s exhibits a counter part in high density molecular gas 
clumps dKoo & Moonl[l9 97b) sharing the same location and ve- 
locity shift. Koo & Moon showed that the shocked gas is dis- 
played in a thin layer in the interface between the SNR shell, 
as delimited by the X-ray image from ROSAT and the un- 
shocked molecular gas. This can be taken as the existence of 
a J-type s hock penetrating the dense gas in a particular region 
of W51B (lKoo & Moonl fl997b). whereas in the location of the 
1720 MHz OH masers the shock s hould be continuous ( C-type). 
Moreover, recent measurements (ICeccarelli et al.ll201 lb showed 
over-ionization of the gas in W51B in certain locations close 
to W51C coinciding with the shocked gas. They conclude this 
excess in ionization implies the existence of an intense flow 
of freshly accelerated cosmic rays (CRs) that, through proton- 
proton collisions, ionize the hydrogen in the adjacent cloud. 
However, ~ 0.2 degrees South-East to the shocked gas region, 
a hard X-ray source CXO J192318. 5+140305 is dete cted. This 
object was first resolve d by ASCA (|Kpo et al.1 12002) and later 
confirmed by Chandra dKoo et al.1 120051) . Its X-ray spectrum, to- 
gether with its morphology, suggests that it is a possible pul- 
sar wind nebula (PWN) associated with the SNR. Therefore, the 
presence of CXO J192318.5+140305 plays a role in the inter- 
pretation of the gamma-ray emission from the W5 1 region. For 
these reasons, W51C represents an interesting case for the study 
of the acceleration of particles to very high energies (VHE) and 
their interaction with the interstellar medium. 

An extended source of gamma rays was first detected by the 
H.E.S.S. telescopes with an integ ral flux above 1 Te V of about 
3% that of the Crab Nebula dFiasson et all l2009h . However, 
the presented morphological and spectral information was not 
enough to attribute the origin of the emission to any particular 
object in the field of view. Also, the Large Area Telescope (LAT) 
on board the Fermi satellite detected an extended source be- 
tween 200 Me V and 50 GeV coincident with the H.E.S.S. source 
dAbdo et al.l l2009al) . Moreover, the reanalysis of the archival 
MILAGRO data after the release of the first Fermi catalog re- 
vealed a 3.4<x excess with medi an energy of 10 Te V coinci- 
dent with the FermifLAT source dAbdo et al.l l2009b). At radio 
wavelengths, synchrotron radiation on ambient magnetic field 
explains the emission detected from W51C. At higher energies, 
there are several processes that yield emission of gamma rays: 
inverse Compton scattering of electrons on seed photons (cosmic 
microwave background, starlight), non-thermal bremsstrahlung 
of electrons on charged target, and decay of neutral pions created 
in flight from a pro ton-nucleon collision. The modeling done by 
( Abdo et al. 2009a) of the spectral energy distribution (SED) of 
W51C disfavors leptonic models and suggests a hadronic ori- 
gin for the emission. For the hadronic channel, two main (non- 
exclusive) mechanisms are to be considered: molecular cloud il- 
lumination by cosmic rays that escap ed the accelerating shock 
dGabici et al.l l2009: Ohir a et al.ll201 ll) or emissi on from clouds 
that are being overtaken by the SNR blast wave dUchivama et al.l 
l2010t iFang & Zhangll2010l) . It is well known that a 10% of the 



energy released by the supernova explosions in the Galaxy can 
account for the energy budget of the CR spectrum up to ener- 
gies close to the knee (~ 10 15 eV). Nevertheless, the evidence 
that SNRs can accelerate particles up to such high energies is 
still missing. Since W51C is one of the most luminous Galactic 
sources at ferm//LAT energies, observation of gamma rays up to 
several TeV would have serious implications regarding the SNR 
contribution to the Galactic CRs: such an observation would 
show that SNRs are not only capable to provide a sufficient flux, 
but could also shed light on the question of the maximum energy 
of CR's still achievable in such a medium age SNR. 

However, the object from which the gamma rays originate 
has not yet been identified within the W5 1 field, and the gamma- 
ray spectrum has so far been precisely measured only up to some 
tens of GeV. In what follows, we report observations with the 
MAGIC telescopes, which will help to address some of the re- 
maining questions on the gamma-ray source in the W5 1 region, 
both regarding its precise location and the physical processes 
needed to explain the observations. In|2]we describe the observa- 
tions that we performed; in[3]we show the observed morphology 
and spectral properties; and, finally, in|4]we apply a theoretical 
framework that can explain the detected gamma-ray emission. 

2. Observations 

MAGIC consists of two 17 m diameter Imaging Atmospheric 
Cherenkov Telescopes (IACT) located at the Roque de los 
Muchachos observatory, on La Palma island, Spain (28°46'N, 
17°53' W), at the height of 2200m a.s.l. The stereo observa- 
tions provide a sensitivitjfl o f 0.8% of the Cr ab Nebula flux at 
energies > 300 GeV, see lAleksic et al.1 d2012l) . MAGIC has the 
lowest trigger threshold of all operating IACTs, enabling it to 
observe gamma rays between 50 GeV and several tens of TeV. 
MAGIC observed W51 in 2010 and 2011. In the first period 
of observations between May 17 and August 19 2010 about 31 
hours effective time remained after quality cuts. Between May 
3 and June 13 2011 additional 22 hours effective time of good 
quality data were taken, resulting in a total amount of 53 h ef- 
fective dark time and covering a zenith angle range from 14 to 
35 degrees. The observations were carried out in the so-called 
wobble mode around the center of the ferm//LAT source W51C 
(RA = 19.385 h, DEC = 14.19°). All data were taken in stereo- 
scopic mode, recording only events which triggered both tele- 
scopes. To minimize systematic effects in the exposure and to 
optimize the coverage for an unknown extension of the emission 
a total of six pointing positions (n po i n t = 6), were used. In all 
pointing positions the wobble distance (offset from the central 
position) was 0.4°, as it is regularly done in MAGIC observa- 
tions. 

The analysis of t he data was performe d using the MARS 
analysis framework dMoraleio et al.l 120091) in cluding the lat- 
est st andard routines for stereoscopic analysis (Lombard i et al.l 
After calibrating the signal and cleaning the images of the 
two telescopes individually, the two images of each stereo event 
are combined. The arrival direction is determined from the com- 
bination of the individual telescope information. To suppress the 
background, a global variable dubbed hadro nness is determine d 
by using the so-called random forest method d Albert et al.l2 008). 
The energy of individual events is estimated using look-up ta- 
bles generated from gamma-ray Monte-Carlo events. For a de- 

1 Sensitivity is defined here as the minimal integral flux to reach 5 cr 
excess in 50 h of observations, assuming a spectral index like that of the 
Crab Nebula. 
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taile d description of the complete analysis chain described above 
see lAleksic etaU (120121) . The gamma-ray signal is estimated by 
comparing the spatial distribution of gamma-like events around 
the assumed source position (ON region) with respect to those 
recorded in signal-free (OFF) regions. The total signal of the 
source is evaluated using a cut on the squared angular distance 
between reconstructed gamma-ray direction and source position 
of 8 2 < 0.07. For each pointing position, the ON sample is com- 
pared to an OFF sample obtained from the combination of the 
"point - 1 OFF regions observed at the same focal plane coor- 
dinates but from the complementary pointing positions. Four of 
the pointing positions have an observation time of the order of 
~ 12 hours each. Therefore, three background samples per point- 
ing can be averaged. The remaining two positions have an ob- 
servation time around 2 hours each and, in this case, the back- 
ground was estimated from one sample only. This method en- 
sures a maximum usage of symmetrical OFF positions without 
introducing big scaling factors due to differences in the obser- 
vation time. The significance of the excess is determined from 
the combined 6 2 dis tribution of all in dividual pointing positions 
using equation 17 in lLi & Mad 1983b (Li&Ma significance here- 
after). 

3. Results 

3.1. Detection 

FigureQ]shows the relative flux mapQ above an energy threshold 
of 150 GeV around the center of the observations. The angular 
resolution of MAGIC for this analysis is 0.085°defined a s one 
sigma of a Gaussian distribution, see lAleksic et al.l d2012l) for 
details. The map was smeared with a two-dimensional Gaussian 
kernel with a sigma equivalent to that of angular resolutiotfl 
Contours represent isocurves of test statistics (TS) evaluated 
from the excess of gamma-like events over a background model. 
This test statistic is Li&Ma significance, applied on a smoothed 
and modeled background estimation. Its null hypothesis distri- 
bution mostly resembles a Gaussian function, but in general can 
have a somewhat different shape or width. The signal region 
is defined within 0.265°radius around the Fermi/LAT position. 
This radius is selected in order to include the emission observed 
in the relative flux map. We compute an excess of 1371.7+ 122.5 
events inside the signal region, yielding a statistical significance 
of 1 1.4 standard deviations. The centroid of the emission (black 
dot in Fig.Q] statistical errors are represented by the ellipse) has 
been derived by fitting a 2 dimensional Gaussian function to the 
map, prior to the smearing. As the centroid we find: 

RA = 19.382 ± 0.001 h DEC = 14.191 ± 0.015 

This deviates by 0.04°from the position reported by Fermi/LAT, 
marked as the center of the sky map (green cross) (see Fig.[TJ. 

To determine the extension of the source we computed the 
distribution of the squared angular distance Q 2 between the ar- 
rival direction of the gamma-like events and the centroid of the 
MAGIC source (see Fig|2]i, both for the integration area repre- 
sented in Fig. [T]and for a combination of signal-free regions from 
where we estimate the background. 

We then fit the difference between ON and OFF 8 2 dis- 
tributions using an exponential function (corresponding to a 

2 Relative flux means excess events over background events. This 
quantity accounts for acceptance differences between different parts of 
the camera 

3 The PSF shown in all skymaps is the sum in quadrature of the in- 
strumental angular resolution and the applied smearing. 




Fig. 1. Relative flux (excess/background) map above 150 GeV 
around W51. Overlaid are contour levels from test statistics 
starting at 3 and increasing by one per contour. The map was 
smoothed with a Gaussian kernel of 0.085°. The green cross rep- 
resents the center of the observations, while the green dashed 
circle represents the integration area. The black dot is the deter- 
mined position of the centroid with the statistical uncertainties 
shown by the surrounding black ellipse. The reg ion of shocked 
atomic and molecular gas dKoo & Moonl 1997 b a) is represented 
by the red dashed ellipse. The blue diamond shows the position 
of the possible PWN CXO J192318.5+140305. In the left lower 
corner the gaussian sigma of a point-like source (PSF) after the 
applied smearing is shown. 



Gaussian-shaped source). For illustration, the shape of a point 
source with the same excess was calculated from Monte-Carlo 
simulations and is shown as comparison (red curve) to the fit to 
the data (blue curve). After correcting for the angular resolution 
(0.085 degrees > 150 GeV) of the instrument the intrinsic exten- 
sion of the source is determined to be: 0.12 + 0.02 stat + 0.02 syst 
degrees. 

3.2. Spectrum 

We extracted the energy spectrum of the gamma-ray emis- 
sion. The effective area was estimated using a Monte-Carlo 
data set with photons simulated uniformly on a ring of 0.15 
to 0.55°distance to the camera center. This accounts for varia- 
tions of the acceptance across the area of the source. The ef- 
fect of using this ring Monte-Carlo compared to standard point- 
like ones turns out to lie well within the statistical uncertainties. 
The spectrum needs to be unfolded in order to take into account 
the fi nite energy resolu tion and the energy bias of the instru- 
ment dAlbert et al.l l2007). The spectrum shown in Fig.[3]starts at 
75 GeV and is well described (^f 2 /NDF = 5.26/6) by a simple 
power law of the form: 
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Fig. 2. 9 2 distribution of the excess events towards the centroid of 
the emission determined from figure [1] showing a clear and ex- 
tended signal. The excess has been fitted by an exponential (blue 
curve) to determine the extension. For comparison the shape of a 
point-like source with the same excess determined from Monte- 
Carlo simulations is shown (red curve). The energy threshold of 
this analysis is 150GeV. 
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Fig. 3. Differential energy spectrum of W51 obtained by 
MAGIC. The red points represent the differential flux points af- 
ter unfolding. The red line represents a power law fit to the data. 
The error bars represent the statistical errors. For comparison, 
the dotte d line represents the spectrum of the Crab Nebula as 
shown in Alek sic et all d2012l) . 



with a photon index of F = 2.58 + 0.07 s t a t ± 0.22 sys t, and 
a normalization constant at 1 TeV of A^o = (9.7 ± 1.0 s tat) x 
10~ 13 cm~ 2 s -1 TeV _1 . This is the first time that the differential 
energy spectrum at VHE is published. The energy threshold 
of MAGIC allows us to almost con nect the spectrum to the 
Fermi/LAT points dAbdo et alj|2009al) . The systematic error on 
the flux normalization is 15%, which includes the systematic un- 
certainties of the effective area (11%) and the background cal- 
culation. In addition, the systematic uncertainty in the energy 
scale is estimated to be 17 % at low (~ 100 GeV) and 15 % at 
medium (~ 250 GeV) energies. The integrated flux above 1 TeV 
is equivalent to ~ 3% of the flux of the Crab Nebula above the 
same energy, and therefore agrees wit h the previous flux mea- 
surement by the H.E.S.S collaboration (Fias son et al.l2 009). The 
spectral index measured by MAGIC a grees well with the o ne 
measured by Fermi/LAI above 10 GeV (Paneque et al. 201 1) of 
T = 2.50 + 0.1 8 s t a t. The emission from W51 can be described by 
a single power law between 10 GeV and 5.5 TeV. 

3.3. Detailed morphology 

MAGIC reaches its best sensitivity in the energy range from 
-300 to -1000 GeV. At energies of 300 GeV the angular reso- 
lution of MAGIC is 0.075 °and it improves until reaching the 
saturation value of 0.054 °at energies above 1 TeV. We investi- 
gate sky maps in two energy ranges. The first map covers the 
estimated energy range from 300 to 1000 GeV, and the second 
the energies above 1000 GeV. Both maps were smeared with a 
Gaussian kernel of a width equal to the angular resolution of the 
instrument in each energy range. 

In Fig. |4](top panels) the relative flux map between 300 and 
1000 GeV is shown. The overall shape of the emission appears 
to be elongated showing a tail towards the lower left. The max- 
imum of the emission coincides with the shocked-gas region, 
represented by the red dashed circle, where the lack of molecu- 
lar material at the systemic velocity is clear (top left panel). The 
determined centroid and extension agree within statistical errors 
with those found above 150 GeV. 



Above 1000 GeV (Fig. [4] bottom panels) the centroid and the 
extension of the emission are in agreement with those obtained 
at lower energies. The South-Eastern tail of the source, evident 
in the 300 to 1000 GeV map, becomes a prominent feature 
coincident with the possible PWN CXO J192318.5+140305 at 
energies above 1 TeV. However, the main part of the emission is 
still coincident with the shocked gas region. 

While the centroid of the emission is consistent with the po- 
sition of the shocked gas, we see a tail towards the PWN can- 
didate. We note that, in any case, the VHE emission does not 
strictly follow the SNR shell (as seen from the 2 1 cm continuum 
emission represented by green contours in the right panels), nor 
does it follow the molecular gas with the velocity expected due 
to Galactic rotation, as traced by the 13 CO (green contours, left 
panels). The tail seen towards the PWN rises the question of a 
possible substructure in the emission. 

3.3.1. Projections 

In order to investigate the source for underlying structures, we 
project the unsmeared excess distribution of the source along a 
line. The line is 2°long divided in 40 bins with 0.05°width. The 
orientation of the line is defined by the position of the PWN can- 
didate and the centroid of the shocked clouds identified by Koo 
& Moon (RA = 19.380 h, DEC = 14.19°). Events within a dis- 
tance of 2 gaussian sigma of the instrumental PSF to the line 
were projected. Since the angular resolution is energy depen- 
dent, the width of the projected rectangle is 0.3°and 0.216°for 
the energy ranges from 300 to 1000 GeV and above 1000 GeV, 
respectively. OFF events were estimated from the background 
model. The number of projected excess events is not the same 
as in the spectral calculation, where we used a circular region of 
0.265°radius around the center of the observations. Therefore, 
the projected excess does not allow for direct determination of 
the fluxes from specific regions of the map. The projection has 
been carried out in both energy ranges independently on the un- 
smeared excess distribution and is shown in Fig. [5] 
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Fig. 4. Relative flux maps: From 300 GeV to 1000 GeV (top) and > 1000 GeV (bottom). On the left hand side the MAGIC data are 
combined with the I3 CO (J=l-0) intensity maps from the Galactic Ring Survey (see http://www.bu.edu/galacticring/newjndex.html) 
integrated between 63 and 72kms~ 1 show n as green countours . On the right hand side the green contours represent the 21cm 
radio continuum emission is shown from (IK00& Moon|[T997ah . In all maps the blue diamond represents the p osition of CXO 
J192318. 5+140305 and the black cross the position of the OH maser emission (|Koo et alj|2005t iGreen et al.l l 1997b . The red dashed 
ellipse represents the region of shocked atomic and molecular gas dKoo & Moonlll997bllah . The 3 counts contour above 1 GeV 
determined by FermifLAI is displayed by the pink contour. In each picture the gaussian sigma of a point-like source (PSF) after the 
applied smearing is shown. The color scale (blue to red) represents the relative flux as measured with MAGIC. In addition the TS 
contours (cyan) are shown starting at 3 and increasing by one per contour. 



We fit the projection alternatively using one and two 
Gaussian functions. ^ 2 /d.o.f. values are 28/17 (one Gaussian) 
and 18/14 (two Gaussians) for the medium-energy range and 
16/17 (one Gaussian) versus 12/14 (two Gaussians) for the high- 
energy events. The data are very well described with the two 
Gaussian functions, where the centroid of the individual func- 
tions coincides within statistical errors with the position of the 
shocked gas and the PWN. The tail-like feature towards the pos- 
sible PWN is more peaked in the energy range above 1000 GeV. 



The statistics are not sufficient to clearly discriminate be- 
tween an extended source of Gaussian excess, an extended 
source of a more complicated shape, or two individual sources. 
However, the fact that there is no region of dense gas close to 
the PWN makes it difficult to explain the enhancement of TeV 
emission in this area under the assumption of uniform CR den- 
sity. A possible scenario of two emission regions could manifest 
in different spectral behaviours. 
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Fig. 5. Projection of the excess inside the marked box in both differential sky maps: 3 00 GeV to lOOOGeV ( top) and above 1000 GeV 
(bottom) along the line connecting the PWN and the shocked-gas region described in lKoo & Moonl(fl9 97a b). The projection is done 
with the unsmeared distribution. The excess is fitted with one (black) and two (red) Gaussian curves. The positions of the shocked 
gas and the PWN are marked with red arrows. On the right-hand side a sketch of the skymaps in both energy ranges is shown to 
illustrate the projected areas, as well as the position of the cloud and the PWN, respectively. The box has a length of l°and a width 
of 4 gaussian sigma of the instrumental PSF. The sky maps show the smeared excess (for comparison with Fig. HJ with the black 
contour representing the 3 TS contour. 



3.3.2. Energy spectra of individual regions 

To quantify the results obtained from the projections we investi- 
gated in more detail the spectral properties of the detected signal, 
we concentrated on two individual regions within the source and 
analyzed them separately. One was defined to cover the shocked 
cloud region with centroid at RA = 19.380h,DEC = 14.19°; 
this will be called the cloud region. The second one was defined 
by the position of CXO J192318.5+140305 and will be called 
the PWN region. To avoid contamination from the surrounding 
emission, and their possible spread due to the worse angular res- 
olution at lower energies, we use an integration radius of 0.1°. 
We compared the same analysis on data of the Crab Nebula and 
find that such a region contains at least 70% of the excess from a 
point-like source above 300 GeV. For an easier comparison, the 
integration radii were chosen to be the same. The distance be- 
tween the chosen positions is 0.19°. There is an area of overlap 
of 1 .7% compared with the integration area of each region, there- 
fore they can be treated as independent. The combined areas of 
both regions represent about 57% of the area used to determine 
the overall spectrum. 

The small distance between the regions and a very similar 
average distance to the camera center allow us to assume the 



same acceptance of gamma-like events for both regions, at least 
within 5%. 

For each individual region we determined the amount of ex- 
cess events above three different energies, and calculated the 
contribution to the overall emission. The resulting values are 
shown in Table [TJ Excesses used to calculate these ratios show a 
significance of at least 2.9cr. 

Table 1. Number of excess events determined for the PWN- 
region and the c/cwaf-region and their contribution to the overall 
emission. Within the statistical errors we do not detect a sig- 
nificant energy dependence on their contributions to the overall 
excess. 



£[GeV] 


cloud 


PWN 


cloud/all [%] 


PWN/an [%] 


> 300 


200 ± 30 


132 ±25 


30 ±5 


19 + 4 


> 500 


116 ± 17 


79 ± 17 


32 ± 6 


22 ±5 


> 1000 


48 ± 10 


27 ± 10 


43 ± 12 


24 ± 10 



The excess contribution arising from the cloud region is 
about 30% and shows no dependence on energy. We performed 
a spectral analysis of a point source for the cloud region above 
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350 GeV. The emission can be well described by a pure power 
law with a flux normalization constant at 1 TeV of V c i olK } = 
(4.3 + 0.9 sta t) x 10" 13 cm- 2 s _1 TeV _1 . The integrated flux above 
350 GeV is equivalent to 1.2% of the flux of the Crab Nebula. 
The spectral index of the cloud emission is -2.4 + 0.5 s tat and 
agrees within statistical uncertainties with the spectral index of 
the overall emission. 

Assuming a point-like emission, the flux from the PWN re- 
gion above 350 GeV is equivalent to 0.7% of the flux of the 
Crab Nebula representing about 20% of the overall observed 
emission. The emission between 350 GeV and 2 TeV can be 
well described by one single power law with a spectral index 
of -2.5 + 0.6 s tat and a flux normalization at 1 TeV of Npwn = 
(2.3 ± 0.8 sta t) x 10- 13 cm- 2 s -1 TeV -1 . 

The excess contribution of each of the regions shows no 
dependence on energy, suggesting no intrinsic morphological 
changes in the energy ranges investigated here. This is in agree- 
ment with the spectra, with the differential maps, and with the 
projections of the excess distribution. We want note that the 
number of excess events within the PWN region and the cloud 
region (Table[TJ agrees within statistical errors with the projected 
excess (Pig. [5j found within +0.1 degree from the PWN and the 
cloud positions, respectively. By looking at the skymaps (Pig. [5j 
only, the emission around the PWN seems to be more intense 
above 1 TeV. However this can be explained by the worse angu- 
lar resolution at lower energies and by a much higher signal-to- 
noise ratio at the higher energies. 



4. Discussion 

Before modeling the multi-wavelength emission we address 
shortly the possibility that the PWN alone is the source of all 
the emission. Then we qualitatively discuss the possibility of a 
contribution of the PWN to the overall emission and justify the 
approach using a one-zone model (i.e. one homogeneous emis- 
sion region filled with one particle distribution per species) to 
investigate the processes underlying the emission. 

First, we want to recall that we have found no spectral or 
morphological energy dependence. In order to assess whether 
the VHE emission can be originated only by the PWN candi- 
date, we consider the estimated rate of rotationa l energy loss 
E- 1 .5 x 10 36 erg s~ ' estimated bvlKoo et all (120051) with the em- 
pirical relation from lSeward & Wangl dl988). Given the obs erved 
lumino sity of the order of ~ 10 36 erg s _1 reported in lAbdo et al.l 
(2009a), it seems unlikely that the PWN alone is the source of all 
gamma-ray emission since, it would require an extremely high 
efficiency in the conversion from rotational energy into gamma 
rays. 

Second, we consider the possibility of a two-zone model. 
The PWN region can account for the 20% of the gamma-ray 
emission confined in a point-like source, however the brightest 
part is the cloud region. This scenario would require an E con- 
version into gamma rays of the order of 10%, in agreement with 
the generally accepted value. 

With the current statistics and resolution, it cannot be es- 
tablished if there is a spectral difference between the cloud and 
the PWN region, but the contribution of the PWN is in any case 
small. For the reasons above, the simplest approach is to assume 
one overall particle distribution underlying the emission we ob- 
serve. This assumption introduces an error in the flux normaliza- 
tion of about 20% in case part of the emission originates from 
the PWN candidate; this uncertainty lies within the statistical 
and systematic errors of the MAGIC measurement. 



4.1. Model description 

We model the SNR as a sphere homogeneously filled with hy- 
drogen, helium and electrons, with respective average number 
densities «h, «He and n e - . For the relative abundances of helium 
we assume the cosmic abundance ratio «He = 0.1 «h- For the 
electron ratio we assume full ionization of the medium, such 

that n e - = 1.2 nu. The magnetic field B is assumed to be ho- 

i i j I 

mogeneous inside the sphere; Koo et al. (2010) derived an upper 

limit for it of B\\ < 150 fiG, but lBrogan et a l. (2000) measured 

a local magnetic field as high as 1.5-1.9 mG towards the maser 

sites. 

The geometric model of the SNR and molecula r cloud in- 
teraction region as proposed by Koo & Moo^ d!997al) . describes 
a scenario in which the spherical blast wave of the supernova 
explosion interacts with part of th e cylindrical molecula r cloud 
contained inside the SNR volume. Carpenter & Sanders ( 1998) 
estimated the total mass of the molecular cloud t o be m c ioud = 
1.9 x 10 s M . From the radio measurements in iMoon &Koo 
(1 1994 ) the angular extent of the partial radio shell of the SNR 
is known to be 9 « 30' . We see a clear displacement between 
the morphology presented here and the center of th e spherical 
extended SNR as seen in thermal X-ray emission (Ko o et all 
1995b). The maximum of the emission is located at the inter- 
action region of remnant and the molecular cloud. Therefore we 
conclude that the size of the remnant is not physically related 
to the size of the VHE emission region. We adopt the intrisnic 
extension determined in this work to determine the radius of a 
spherical emission z one. Assuming a dist ance to W51C of 5. 5 
kpc, as measured bv lSato etaTI (1201 Oh and lMoises et all d2011l) . 
the radius of the sphere is estimated to be 24 pc. 

The explosion energy of the SNR has been estimated in 
iKoo et all d!995bl) as E SN « 3.6 x 10 51 erg, using both a Sedov 
and an evaporative model to derive the parameters of the SNR. 
We will compare this value with the one obtained from the in- 
tegral of our initial spectra, after we fix the normalisation con- 
stants; we will determine how much of the initial explosion en- 
ergy of the supernova has been converted into particles ( W e , W p ). 
The different parameters of the supernova, of the SNR, and of the 
molecular cloud are summarized in Table|2] 



Table 2. Parameters of the W51C supernova, supernova remnant 
and molecular cloud. 



Parameter 


Value 


Reference 


age 


~ 30 000 yr 


Koo et al. (1995b) 




~ 3.6 x 10 51 erg 


Koo et al. ( 1995W 


d 


5.5 kpc 


Sato et al. (2010) 






Moises etal. (2011) 


9 (radio) 


S3 30' 


Moon & Koo (1994) 




< 150 fiG 


Koo et al. (2010) 


B (at masers) 


1.5-1.9 mG 


Broean et al. (2000) 


a T 


S3 -0.26 


Moon & Koo (1994) 


m cloud 


1.9 x 10 5 M G 


Carpenter & Sanders (1998) 



We model the spectral energy distribution folding input spec- 
tra of accelerated particles with cross sections of processes 
yielding photons; this includes synchrotron radiation, inverse 
Com pton scattering (IC), non-thermal bremsstrahlung and n° de- 
cay dBlumenthal & Gouldll970HBaring et all 1 9991; iKelner et all 
2006). 

For IC, we consider three seed photon fields: the cosmic 
microwave background (IcTcmb = 2.3 x 10~ 4 eV, mcmb = 
0.26eVctrT 3 ), infrared (kT 1R = 3 x 10~ 3 eV, u 1R = 0.90eVcm" 3 ) 
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and optical (kTopj = 0.25 eV, «opt = 0.84 eV cm" 3 ), with 
temperatures and energy densities for the inf rared and optical 
components adopted from lAbdo et al.l (l2009al) . Bremsstrahlung 
is computed on a target of electrons and ions. For the n° pro- 
duction cross section, we use the parametrization of Kelner & 
Ahar onian (200 6) with a constant nuclear enhancement factor of 
1.85 (lMorill2009l) . 

The multi-wavelength data considered he re include radio 
continuum measurements (Moon & Koo 1994), high-energy ob- 
servations by the Fermi/LAT dAbdo et all l2009ah and the new 
VHE data taken with MAGIC, present ed in this paper. In cluded 
is also one data point by MILAGRO (Ab do et al.ll2009bl) . Note 
that the lowest ener gy radio data point m ay be affected by free- 
free ab sorption, see Moon &K00I (1 1994 or lCopetti & Sc hmidt 
(fl99l . which we do not consider here. However, this single 
point does no t affect the fitting of t he radio data. The radio mea- 
surements in Moon & Kool (11994 indicate a spectral index of 
a,- « -0.26 (as defined by S v oc v" r ). This can be attributed 
to electrons emitting synchrotron radiation and fixes the initial 
power-law index of the electron spectrum to s « 1.5. We adopt 
this value both for electrons and protons. 

As an upper limit for the non-thermal X-ray emission we 
consider the integrated t hermal X-ray flux of the whole remnant 
as measured by ROSAT (Koo et aDl995bl) converted into a dif- 
ferential flux in the sub-keV range. We use the thermal emis- 
sion observed by Chandra from CXOJ192318.5+140305 as an 
upper limit to the non-thermal emission of the possible PWN. 
The MILAGRO measurement has a significance of 3.4 cr, was 
derived assuming a gamma-spectrum oc E~ 2 (> witho ut a cut-off 
and is g iven at an energy of 35 TeV. For details see Abd o et aT] 
(l2009bl) . 

We consider seperate scenarios in which one of the following 
emission processe dominates over the others, pion decay, inverse 
Compton, or Bremsstrahlung. The models discussed here are ob- 
tained using as equilibrium particle spectra a broken power law 
with an exponential cut-off, both for electrons and protons, of 
the form: 



dN, 



e.p 



dE, 



e.p 



£ 



1 + 



Ebr 



As 



exp 



J cut,e,/7 



(2) 



The spectral index changes here from s to s + As at an energy 
£br with a smooth transition. The exponential cut-off at E cat ^ p 
reflects the roll-off of the particle spectrum near the maximum 
energy, arising from the acceleration and confinement mecha- 
nism, as well as energy losses. 

The break energy £\,r is fixed from the Fermi/LAT data, 
while the new MAGIC data allow us to fix the spectral break 
As. A spectral break in the particle spectrum at these en- 
ergies is traditionally thought to be inconsistent with both 
standard or non-linear diffus ive shock acceleration theory, see 
Malk ov & O'C Drurvl (l200lh and references therein. However, 
Malk ov et al.l (1201 ll) have recently proposed a mechanism which 
can also explain a spectral break in the cosmic ray spectrum of 
As — 1 by strong ion-neutral collisions in the surroundings of a 
SNR, leading to a weakening in the confinement of the acceler- 
ated particles. The spectral break that we have derived here is As 
= 1.2, not far off this prediction, giving a hint that this mecha- 
nism might be responsible for the observed break. Note also that 
other authors have proposed scenarios in which the CR spec- 
trum, and consequently the gamma-ray spectrum, can show one 
or more spectral breaks, for example due t o finite-size accelera- 
tion or emission region (Ohira et al. 201 1) or energy dependent 



diffusion of run-away CRs from the remnant dGabici et al.l t2009: 
lAharonian & AtovarJl 19961) . 

The luminosity of W51C in the energy range 0.25 GeV - 
5.0 TeV, which is roughly the energy range of the Fermi and 

assuming a distance of 



MAGIC data, is L y 



1 x lO^ergs" 1 , 



5.5 kpc, which is one of the highest compared with other SNRs. 



4.2. Adjustment of model parameters 

First we consider the case where the emission is dominated 
by leptonic emission mechanism s . We fi nd the same problems 
already reported by lAbdo et aD d2009al) . namely that it can- 
not reproduce the radio and gamma-ray data simultaneously. 
Furthermore these models need an unusually high electron to 
proton ratio of the order of one. 

When we model the emission with pion decay as the dom- 
inant process, both radio and gamma-ray emission can be rea- 
sonably reproduced, as shown in Fig. [6] A hadronic scenario 
is particularly interesting, as the shock-cloud interaction natu- 
rally favors a CR-matter interaction mechanism. Moreover, the 
parameters used in this model, see Table [3] are a reasonable de- 
scription of the interstellar medium around W5 1 . 



Table 3. Parameters used in the modeling of the multi- 
wavelength spectral energy distribution for the hadronic sce- 
nario. The power-law index before the break is s = 1.5 for both 
protons and electrons. Eq — 10 GeV. The total kinetic energy of 
the particles was integrated for E\^ n > 100 MeV both for elec- 
trons and protons. 



Parameter 



Value 



K e /K p 


1/80 


As 


1.2 


E bl [GeV] 


10 


E M , e [TeV] 


0.1 


E cutp [TeV] 


120 


B\mG] 


53 


n [cm -3 ] 


10.0 


W e [10 50 erg] 


0.069 


W„ [10 50 erg] 


5.8 



Compared to th e hadronic model suggested in the work of 
lAbdo et alJ d2009a). the main difference is the index of the part- 
cile distribution after the break. The spectrum after the break is 
more precisely determined by the data presented here. The in- 
dex we obtain is harder, allowing for the explanation of all the 
gamma-ray data up to the end of the MAGIC spectrum. 

A detailed view of the high energy and VHE region is shown 
in Fig. [7] It shows that the index above the break is clearly de- 
termined by the data presente d here. In addition, the hadronic 
model by lAbdo et ail d2009al) is displayed. In addition to the 
good aggrement between the model and the data, the plot shows 
that the results presented here clearly improve the determination 
of the underlying particle distribution. In this scenario a cut-off 
energy of E cutp > 100 TeV is needed to fit the MAGIC data, 
indicating the existence of protons at least to this energy. 

The precise cut-off energy of the electron spectrum £ culf is 
not well constrained, since the synchrotron peak is not resolved. 
Therefore, the energy E cut<e used here is only a lower limit, as 
enforced with the radio data. However, a 1 TeV electron in a 
magnetic field of 50 /iG has a lifetime of about 4700 years, de- 
termined by synchrotron losses. This value is much lower than 
the age of the remnant, suggesting that for such high energies 
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Fig. 6. Model of the multi- wavelength SED in the hadronic - 
dominated scenario. The dashes with error bars are 21cm ra- 
dio continuum, circles represent fen«//LAT data, squares are the 
data obtained in this work and the star represent the MILAGRO 
data point. The upper limit in the X-ray regime is obtained from 
ROSAT data as discussed in the text. The details of the scenario 
are discussed in the text. 
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Fig. 7. Detailed view of the hadronic model in the high en- 
ergy and VHE regio n. For comparison to the hadronic model by 
Abdo et al. (2009a), shown as double dotted line. This models 
ends at the highest energy shown in that p ublication. The ma- 
jor difference between our model and that of lAbdo et alJ (|2009a) 
is the harder particle spectrum above ~100GeV, which is now 
precisely constrained by the measurements presented here. 



The volume-averaged hydrogen density is obtained as a 
parameter of the fit. From that, we compute the volume fill- 
ing factor /, which is the fraction of the mass of the clumpy 
molecular cloud that is contained inside the SNR interaction 
volume (defined as volume of the emission zone) as / = 
n H v{jn^ om + 0.1mg° m ) /m c i oud * 0.11. Here m c i oud is the total 
mass of the molecular cloud, V is the volume of the radiation 
sphere and m^ om , mS? m are the masses of a hydrogen and helium 
atom, respectively. This would imply that around 1 1 % of the 
mass of the molecular cloud is contained in the emission volume 
and is interacting with the SNR. This value is consistent with the 
filling factors of around 8-20 % for other SNRs interacting with 
molec ular clouds, obtained by other authors dUchivama et alJ 
I2010I) . 

The total amount of kinetic energy in electrons and protons 
is about 16 % of the explosion energy of the supernova. This 
fraction is just slightly higher than the value normally assumed, 
of around 10 %, of the explosion energy con verted into C Rs to 
maintain the observed flux of Galactic CRs (Hillas 2005). The 
proton to electron ratio is not far from value o bserved at earth of 
K p IK e « 50, see for example [Simpsonl d!983l) . 

Since the hadronic gamma-ray emission is proportional to 
the product of the kinetic energy in protons and the density of the 
medium, this parameters are striclty correlated. Assuming that 
the complete mass of the molecular cloud acts as target material 
(f=l), this would imply a density of n=100 cm 3 . Therefore the 
lower limit of the energy in relativistic protons is about 1.6% 
of the explosion energy of the supernova. We note that such 
a scenario would need either a higher magnetic field (B ~ 150 
fiG) or a much lower electron to proton ratio (K e /K p ~ 1/800) 
to still reproduce the broadband emission. In addition, the mor- 
phology presented in this work shows that only a fraction of the 
molecular cloud is emitting VHE gamma emission (see Fig. |4j». 
Therefore we conclude that the amount of kinetic energy in pro- 
tons is clearly above this lower limit and in the order of 10-20%. 

In the scenario investigated here all of the gamma-ray emis- 
sion was attributed to tt° decay. It was not possible to model the 
broad-band emission with a purely leptonic scenario. The radio 
data could not be fitted and the model parameters were not phys- 
ically reasonable (too low density tin, too high energy content 
W e in electrons, too low magnetic field B). However, that could 
also point to problems in the modeling, especially to oversimpli- 
fications concerning the homogeneity of the medium and of the 
magnetic field. 

We conclude that the Fermi/LXT data and the MAGIC data 
can be explained in terms of hadronic interactions of high-energy 
protons with the molecular cloud and subsequent decay of neu- 
tral pions. With the current data it is not possible to decide what 
process causes the hint of emission observed by MILAGRO 
which, if confirmed at this flux level, would require the intro- 
duction of an additional component at the highest energies. 



the electron spectrum should develop a break, with the conse- 
quent spectral steepening at higher energies. Assuming constant 
electron injection over time, the electron spectrum steepens at 
100 GeV by a factor l/E. This yields a very similar gamma- 
ray emission as in the hadronic model presented here, even for a 
higher value of E cuUe . 

4.3. Physical outcome of the models 

We discuss what general conclusions can be drawn from the 
model which fit the data: the hadronic scenario. 



4.4. Discussion on the acceleration process 

Following the result of the modeling we assume the observed 
gamma-ray emission to be of hadronic origin. As mentioned in 
SectionQ] there are two main possible scenarios: a cloud illumi- 
nated by runaway CRs or acceleration of CRs in the shock wave 
propagating through the cloud. 

In the first case, CRs escaping the SNR will homogeneously 
fill a sphere with a radius R c { ~ y4Dt where D is the dif- 
fusio n coefficient and t is the time since particles are diffus- 
ing dGabici etaU r2010). For a distance of 5.5 kpc and lOTeV 



9 



J. Aleksic et al.: Morphological and spectral properties of the W51 region measured with the MAGIC telescopes 



protons, responsible for gamma-ray emission of 1 TeV, the ra- 
dius of this sphere would be about 350 pc, assuming the average 
Galactic CR diffusion coefficient at 10 TeV to be ~ 3xl0 29 cm 2 s. 
Here we assumed that the high-energy particles escape the SNR 
early enough such that the diffusion time can be approximated 
to be the age of the SNR. The distance between the maximum of 
the emission measured by MAGIC above 1 TeV and the assumed 
center of the SNR (RA=19.384 h, DEC=14.11°) is about 8 pc. 
The distance to other parts of the SNR/cloud complex W51C/B 
is of similar order. This implies that the complete cloud should 
be uniformly illuminated by CRs. As can be seen in Fig. [4] we 
do not detect the complete W51B/C complex at energies above 
1 TeV (lower right skymap): parts of the outer regions, both on 
the side towards the SNR and on the opposite side, do not emit 
gamma radiation. In the scenario of runaway CRs, we would 
expect diffusion from the SNR to W51A (northern region in the 
21cm emission); no significant emission from W51A is detected. 
However, the distance of the regions A and B is measured with 
an error of the order of hundreds of parsecs, which means that 
the relative distance between the two could be high enough to 
explain the lack of diffusion from one to the other. The scenario 
of runaway CRs also can not explain the incomplete illumination 
of W51B/C, especially towards the outer regions. 

Concerning the acceleration of CRs in the shocked cloud 
scenario, the gamma radiation should be originated very close to 
the acceleration site of the radiating particles due to the high den- 
sity of the surrounding medium. This is in agreement with the 
morphology des cribed in this work. Th e unusually high ioniza- 
tion reported by ICeccarelli et al.l d201 ll) close to the maximum 
VHE emission region indicates the presence of freshly accel- 
erated low-energy protons. The missing emission towards the 
edges of the cloud could be explained with a lower diffusion co- 
efficient in the shocked cloud region, or with a shielding effect, 
either of which is possible in a surrounding medium of high den- 
sity. 

Both the morphology at TeV energies and the measured high 
ionization are hints for an ongoing acceleration. This suggests 
that the particle distribution, whose gamma emission we ob- 
serve, may represent the source spectrum of cosmic rays cur- 
rently being produced in W5 1 . However, the differentiation be- 
tween ongoing acceleration of particles in the shocked region or 
reacceler ation of already exist ing CRs, like in the crushed cloud 
scenario ( Uchivam a et alj|2010h . in the same region is not obvi- 
ous and is not addressed in this work. 



5. Conclusions 

MAGIC has performed a deep observation of a complex Galactic 
field containing the star-forming regions W51A and W51B, the 
SNR W51C and the possible PWN CXO J192318.5+140305. 
As a result of this observation, emission of gamma rays above 
150 GeV has been detected with 11 cr statistical significance. 
The spectrum of this emission has been measured between 
75 GeV and 4 TeV. Spectral points are well fitted with a power 
law with a photon index of 2.6, compatible with the Fermi/LAT 
measurement between 2 and 40 GeV. The spectrum measured by 
MAGIC allows for the first time a precise determination of the 
spectral slope of the underlying particle distribution above the 
spectral break measured at around a few GeV by Fermi/LAT. 

The MAGIC source spatially coincides with those previously 
reported by H.E.S.S. and Fermi/LAT. We are able to restrict the 
emission region to the zone where W51C interacts with W51B 
and, in particular, to the region where shocked gas is observed. 



This clearly pinpoints the origin of the emission to the interac- 
tion between the remnant and the molecular cloud. 

Non-thermal X-ray emission which could help to trace 
the relativistic electron distribution was found only from a 
compact region around the position of the possible PWN 
CXOJ192318.5+140305 dKoo et alJl2005h . The MAGIC source 
exhibits a morphological feature extending towards CXO 
J192318. 5+140305, more prominent in the image at higher en- 
ergies. 

The projection of the gamma-like events on the line connect- 
ing the putative PWN and the centroid of the shocked clouds 
shows a hint of an underlying distribution that may be described 
as the sum of two Gaussian functions. However, the existence 
of two independent, resolved sources cannot be statistically es- 
tablished. We thus investigate the contribution to the total excess 
of two regions of 0.1°radius centered on the cloud region and 
the PWN region. We find that they contribute about 30% and 
20% of the total emission, respectively, and the contribution is 
not energy dependent within the uncertainties. Spectra of the in- 
dividual regions above 350 GeV could be obtained, but do not 
allow for detailed conclusions due to the weak individual fluxes. 
Given the small possible contribution of the PWN candidate in 
the energies investigated in this work, it is very unlikely that the 
main conclusion drawn here will be significantly affected even if 
the PWN contribution can be established. 

MAGIC observations determine the VHE spectral energy 
distribution of W5 1 over more than one order of magnitude in 
energy. We have produced a physically plausible model of the 
emission of the SNR by considering a spherical geometry and 
uniform distribution of the ambient material. We note that this 
system is clearly anisotropic (as seen in the multi-wavelength 
data), and more detailed modeling may achieve a better descrip- 
tion of the source. We find that the VHE emission from W51C 
cannot be explained by any of the considered leptonic models. 
The emission is best described when neutral pion decay is the 
dominant gamma-ray production mechanism. In the proposed 
model, the SNR has converted about 16% of the explosion en- 
ergy into kinetic energy for proton acceleration and the emission 
zone engulfs a 10% of a molecular cloud of 10 5 solar masses, 
which provides the target material. In this scenario, protons are 
required to reach at least an energy of the order of 100 TeV to 
produce the observed emission. 

The morphology of the source cannot be explained by CRs 
diffusing from the SNR to the cloud. It can instead be qualita- 
tively explained with VHE gamma-ray emission being produced 
at the acceleration site of CRs. This involves ongoing acceler- 
ation of CRs or re-acceleration of already existing CRs at the 
shocked cloud region. Given the high luminosity of this source 
and its plausible hadronic origin, we conclude that W51C is a 
prime candidate cosmic ray source in the Galaxy. 

Finally, we want to give a short outlook and address the a few 
issues connected with W51C. The detection of neutrinos from 
this source would be the final proof about the hadr onic nature 
of the emission. But, according to the calculations bv lYuan et al.l 
(1201 lh . the chances for detection are low. However, also an ex- 
tension of the high-energy gamma emission towards low er en- 
ergies, as performed for example in Giulia ni et al.l d201 lh . may 
also provide more clues to the nature of particle acceleration in 
this region. To reveal the morphology and the possible emis- 
sion of the PWN, more data at energies above 1 TeV are neces- 
sary. Extension of the spectrum towards higher energies would 
constrain the maximum achievable energy in the system and 
might shed light on the meaning of the MILAGRO measure- 
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ment, which cannot be accommodated in the theoretical frame- 
work proposed here. 
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